
#include "../common/modularity_mat.h"
#include "../common/error.h"
#include "../common/mem.h"
#include "../common/matrix_op.h"

#include <float.h>

ModularityMatrix* make_modularity_mat(sparse_matrix_arr* matrix, int* indices,
									  int indices_count)
{
	ModularityMatrix* Bg = mem_alloc_and_zero(1, sizeof(ModularityMatrix));
	sparse_matrix_arr* adj_sub_matrix;
	int i,j;
	int n = matrix->n;
	elem* Fg;
	elem M;
	elem* Kg;
	elem norm1 = -DBL_MAX;
	elem sum;
	int last_summed_row;

	Kg = mem_alloc(indices_count, sizeof(elem));
	Fg = mem_alloc(indices_count, sizeof(elem));

	validate_allocation(Kg, "Could not allocate internal structures!");
	validate_allocation(Fg, "Could not allocate internal structures!");

	Bg->Fg = Fg;
	Bg->Kg = Kg;

	adj_sub_matrix = sparse_matrix__make_sub_matrix(matrix, indices,
		indices_count);
	Bg->Ag = adj_sub_matrix;

	M = 0;

	last_summed_row = -1;
	for (i = 0; i < indices_count; i++)
	{
		/*  TODO: Can be optimized (can be done once for the whole run) */
		Kg[i] = sparse_matrix__sum_row(matrix, indices[i]);

		/*  Get maximal "Kg" value (which is the norm1 value) 
		if (Kg[i] > max_Kg_value)
			max_Kg_value = Kg[i]; */

		M += Kg[i];

		/*  TODO: Can be optimized (calculate "M" once for the original matrix) */
		/*  Sum into M the excluded rows */
		for (j = last_summed_row + 1; j < indices[i]; j++)
			M += sparse_matrix__sum_row(matrix, j);

		last_summed_row = indices[i];
	}
	
	/*  Sum into M the excluded rows */
	for (j = last_summed_row + 1; j < n; j++)
		M += sparse_matrix__sum_row(matrix, j);

	Bg->M = M;
	/* Bg->Norm1 = max_Kg_value; */

	/*  Calculate Fg for the generalized modularity matrix */
	for (i = 0; i < indices_count; i++)
	{
		Fg[i] = 0;

		for (j = 0; j < indices_count; j++)
			Fg[i] += modularity_mat_get_cell(Bg, i, j, 0);
	}
	
	for (i = 0 ; i < indices_count ; i++)
	{
		sum = 0;
		for (j = 0 ; j < indices_count ; j++)
		{
			sum += ABS( modularity_mat_get_cell(Bg, i, j, 1) );
		}
		norm1 = (sum > norm1) ? sum : norm1;
	}
	Bg->Norm1 = norm1;

	return Bg;
}


elem modularity_mat_get_cell(const ModularityMatrix* B, int row, int col,
							 int generalized)
{
	elem result = sparse_matrix__get_cell(B->Ag, row, col);

	result -= B->Kg[row] * B->Kg[col] / B->M;

	if (row == col && generalized)
		result -= B->Fg[row];

	return result;
}

void free_modularity_mat(ModularityMatrix* B)
{
	free_sparse_matrix_arr(B->Ag);
	mem_free(B->Fg);
	mem_free(B->Kg);
	mem_free(B);
}


void modularity_mat_mlpl_by_vector(const ModularityMatrix* Bg, Vector* x,
								   Vector* result, int is_generalized_matrix)
{
	elem scalar;
	Vector* helper;
	int n = Bg->Ag->n;
	Vector dummy;

	helper = allocate_Vector(n);

	/*  result = A[g] * x */
	sparse_matrix__mlpl_matrix_times_vector(Bg->Ag, x->data, result->data);

	/*  helper = (K[g] * x * K[g]) / M */
	dummy.size = n;
	dummy.data = Bg->Kg; /*  Kg */
	scalar = scalar_mlpl(&dummy, x) / Bg->M ; /*  scalar = (K[g] * x) / M */
	scalar_vector_mlpl(scalar, &dummy, helper);    /*  helper = scalar * K[g] */

	/*  result -= helper  < result=A[g]*-((K[g]*x)/m)*K[g] > */
	add_vectors(result, helper, result, SUBTRACT_VECTORS);

	if (is_generalized_matrix)
	{
		/*  helper = D[g] * x  =>  helper(i) = Fg(i) * x(i) */
		dummy.size = n;
		dummy.data = Bg->Fg; /*  f(g) */
		diagonal_matrix_vector_mlpl(&dummy, x, helper);

		/*  result -= helper  < result=A[g]*x-((K[g]*x)/m)*K[g]-D[g]*x > */
		add_vectors(result, helper, result, SUBTRACT_VECTORS);
	}

	free_vector(helper);
}